Projet télédétection / Master 2 ECOMONT
Ce projet est réalisé dans le cadre du Master 2 ECOMONT et porte sur l’étude spatiale des données SMOD (date de fonte de la neige) et NDVI (indice de végétation par différence normalisée) obtenues sur une zone délimitée, située sur le territoire de la communauté de communes de la vallée de Chamonix.
1 Chargement des librairies
library(raster)
library(dplyr)
library(mblm)
library(tiff)
library(rgdal)
library(rgeos)
library(terra)
library(basemaps)
library(mapedit)
library(rasterVis)
library(RColorBrewer)
library(tidyverse)2 Sources des données
2.1 Données NDVI
Institution(s):
Pôle THEIA du CNES (traitement des images
Sentinel-2)
CREA Mont-Blanc (calcul du NDVI)
Auteurs
Pôle THEIA - https://theia.cnes.fr/atdistrib/rocket/#/home
Brad
Carlson (CREA Mont-Blanc)
2.2 Données CCVCMB
Institution(s): Communauté de Communes de la Vallée de Chamonix Mont-Blanc (CCVCMB)
2.3 Données SMOD
Données initiales
CESBIO Toulouse
Pôle THEIA du CNES
Accès : https://theia.cnes.fr/atdistrib/rocket/#/home
Citation : https://www.mdpi.com/2072-4292/12/18/2904/pdf
Contact : simon.gascoin@cesbio.cnes.fr
Traitement des
données :
CREA Mont Blanc
bcarlson@creamontblanc.org
3 Objectifs du projet
Notre projet a consisté à étudier s’il existait une corrélation entre l’évolution du SMOD et l’évolution du verdissement observé (via le calcul du paramètre NDVI), entre 2017 et 2021, sur une zone délimitée de la communauté de communes de la vallée de Chamonix (CCVCMB).
4 Analyse spatiale et statistique des données
4.1 Import et description des couches raster
4.1.1 Import des couches NDVI par année
rastlist.2017.NDVI <- list.files(path = "S2_NDVI_CCVCMB", pattern='2017', all.files=T, full.names=T)
rasters.2017.NDVI <- stack(lapply(rastlist.2017.NDVI, raster))
rastlist.2018.NDVI <- list.files(path = "S2_NDVI_CCVCMB", pattern='2018', all.files=T, full.names=T)
rasters.2018.NDVI <- stack(lapply(rastlist.2018.NDVI, raster))
rastlist.2019.NDVI <- list.files(path = "S2_NDVI_CCVCMB", pattern='2019', all.files=TRUE, full.names=T)
rasters.2019.NDVI <- stack(lapply(rastlist.2019.NDVI, raster))
rastlist.2020.NDVI <- list.files(path = "S2_NDVI_CCVCMB", pattern='2020', all.files=TRUE, full.names=T)
rasters.2020.NDVI <- stack(lapply(rastlist.2020.NDVI, raster))
rastlist.2021.NDVI <- list.files(path = "S2_NDVI_CCVCMB", pattern='2021', all.files=TRUE, full.names=T)
rasters.2021.NDVI <- stack(lapply(rastlist.2021.NDVI, raster)) 4.1.2 Caractéristiques des couches NDVI
## class : RasterStack
## dimensions : 2636, 2348, 6189328, 30 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 324930, 348410, 5077160, 5103520 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs
## names : SENTINEL2A_20200224_NDVI, SENTINEL2A_20200315_NDVI, SENTINEL2A_20200325_NDVI, SENTINEL2A_20200404_NDVI, SENTINEL2A_20200409_NDVI, SENTINEL2A_20200424_NDVI, SENTINEL2A_20200504_NDVI, SENTINEL2A_20200509_NDVI, SENTINEL2A_20200519_NDVI, SENTINEL2A_20200524_NDVI, SENTINEL2A_20200529_NDVI, SENTINEL2A_20200603_NDVI, SENTINEL2A_20200623_NDVI, SENTINEL2A_20200708_NDVI, SENTINEL2A_20200713_NDVI, ...
## min values : -0.9983560, -0.9981877, -0.9985891, -0.9978938, -0.9985730, -0.9980031, -0.9985617, -0.9979624, -0.9981035, -0.9982855, -0.9980518, -0.9969445, -0.9991574, -0.9974279, -0.9978781, ...
## max values : 0.7155070, 0.7039071, 0.7296859, 0.8276377, 0.8230168, 0.8353869, 0.8483894, 0.9216663, 0.8404192, 0.8615274, 0.8427536, 0.8281259, 0.8500417, 0.8184949, 0.8514056, ...
Les rasters des données NDVI issues du satellite Sentinel2 ont une résolution de 10 m par 10 m et sont projetés dans le référentiel SCR WGS84 / UTM zone 32N.
4.1.3 Visualisation d’une couche NDVI (24/02/2020)
Figure 1: Couche NDVI limitée à la CCVCMB le 24/02/2020
Les zones blanches correspondent à des nuages. Les couleurs les plus vertes dessinent la vallée de Chamonix et le gradient des couleurs jaunes à roses refléte le gradient altitudinal des versants qui encadrent la vallée.
4.1.4 Import des couches SMOD par année
SMOD_17 <- raster("Sentinel2_neige/SMOD_CCVCMB_2017.tif")
SMOD_18 <- raster("Sentinel2_neige/SMOD_CCVCMB_2018.tif")
SMOD_19 <- raster("Sentinel2_neige/SMOD_CCVCMB_2019.tif")
SMOD_20 <- raster("Sentinel2_neige/SMOD_CCVCMB_2020.tif")
SMOD_21 <- raster("Sentinel2_neige/SMOD_CCVCMB_2021.tif")4.1.5 Caractéristiques des couches SMOD
## class : RasterLayer
## dimensions : 1318, 1175, 1548650 (nrow, ncol, ncell)
## resolution : 20, 20 (x, y)
## extent : 324920, 348420, 5077160, 5103520 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs
## source : SMOD_CCVCMB_2017.tif
## names : SMOD_CCVCMB_2017
## values : 0, 243 (min, max)
Les rasters des données SMOD issues du satellite Sentinel2 ont une résolution de 20 m par 20 m et sont projetés dans le référentiel SCR WGS84 / UTM zone 32N. Il faut donc redimensionner les couches SMOD pour avoir la même résolution que les couches NDVI (10 m x 10 m).
4.1.6 Changement de la résolution des couches SMOD
SMOD_17 <- resample(x = SMOD_17, y = rasters.2020.NDVI[[1]], method = "bilinear")
SMOD_18 <- resample(x = SMOD_18, y = rasters.2020.NDVI[[1]], method = "bilinear")
SMOD_19 <- resample(x = SMOD_19, y = rasters.2020.NDVI[[1]], method = "bilinear")
SMOD_20 <- resample(x = SMOD_20, y = rasters.2020.NDVI[[1]], method = "bilinear")
SMOD_21 <- resample(x = SMOD_21, y = rasters.2020.NDVI[[1]], method = "bilinear")
SMOD_17## class : RasterLayer
## dimensions : 2636, 2348, 6189328 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 324930, 348410, 5077160, 5103520 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs
## source : memory
## names : SMOD_CCVCMB_2017
## values : -22.25, 287 (min, max)
La résolution des couches SMOD a bien été modifiée.
4.1.7 Visualisation d’une couche SMOD (2017)
Figure 2: Couche SMOD en 2017 autour de la vallée de Chamonix
Des couleurs claires roses/jaunes sont associées à la vallée de Chamonix et donc à des zones pour lesquelles la neige fond précocement. Les couleurs vertes vives sont associées à des zones hautes en altitude pour lesquelles la neige fond tardivement ou jamais.
4.1.8 Application du masque CCVCMB aux couches SMOD
4.1.8.1 Import des limites de la communautés de communes (CCVCMB)
La couche correspondant aux limites de la CCVCMB est importée au préalable pour appliquer son emprise sur les couches SMOD afin d’obtenir des limites similaires aux couches NDVI.
CCVCMB <- shapefile("Communes/CC_Vallee_CHAMONIX.shp")
CCVCMB <- spTransform(CCVCMB, CRS = crs("+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs "))La couche CCVCMB est également reprojetée dans le même référentiel SCR que les couches NDVI.
4.1.8.2 Application du masque CCVCMB
SMOD_17M <- mask(SMOD_17, CCVCMB)
SMOD_18M<- mask(SMOD_18, CCVCMB)
SMOD_19M<- mask(SMOD_19, CCVCMB)
SMOD_20M<- mask(SMOD_20, CCVCMB)
SMOD_21M<- mask(SMOD_21, CCVCMB)4.1.8.3 Visualisation de l’effet du masque CCVCMB
Figure 3: Visualisation de la couche SMOD 2017 après application du masque CCVCMB
4.2 Préparation des tableaux NDVI max et SMOD
L’objectif est désormais de calculer le NDVI max pour chaque année afin de pouvoir calculer la pente liée à l’évolution du NDVI max sur les cinq années d’étude (2017-2021). La pente associée à l’évolution des couches SMOD sur les cinq mêmes années sera également calculée.
4.2.1 Calcul NDVI max par année
Les rasters et les tableaux associés aux NDVI max sont générés.
NDVImax.2017 <- max(rasters.2017.NDVI, na.rm=T)
NDVImax.2017_df <- as.matrix(NDVImax.2017)
NDVImax.2018 <- max(rasters.2018.NDVI, na.rm=T)
NDVImax.2018_df <- as.matrix(NDVImax.2018)
NDVImax.2019 <- max(rasters.2019.NDVI, na.rm=T)
NDVImax.2019_df <- as.matrix(NDVImax.2019)
NDVImax.2020 <- max(rasters.2020.NDVI, na.rm=T)
NDVImax.2020_df <- as.matrix(NDVImax.2020)
NDVImax.2021 <- max(rasters.2021.NDVI, na.rm=T)
NDVImax.2021_df <- as.matrix(NDVImax.2021)4.2.2 Visualisation NDVI max par année
Figure 4: Couches NDVI max par année
Peu de variations sont observées sur les 5 années. Un léger front de verdissement est observé sur le versant est de Chamonix mais il reste compliqué d’identifier des changements facilement observables.
4.2.3 Homogénéisation du nombre de pixels entre les 5 années pour les couches SMOD et NDVI max
Les masques des couches présentant les nombres de pixels différents de “NA” les plus faibles sont appliqués à l’ensemble des couches afin d’obtenir un nombre de pixels similaire pour l’ensemble des couches.
NDVImax.2017<-mask(NDVImax.2017, NDVImax.2020)
NDVImax.2018<-mask(NDVImax.2018, NDVImax.2020)
NDVImax.2019<-mask(NDVImax.2019, NDVImax.2020)
NDVImax.2020<-mask(NDVImax.2020, NDVImax.2020)
NDVImax.2021<-mask(NDVImax.2021, NDVImax.2020)SMOD_17M <- mask(SMOD_17M, SMOD_18M)
SMOD_18M<- mask(SMOD_18M, SMOD_18M)
SMOD_19M<- mask(SMOD_19M, SMOD_18M)
SMOD_20M<- mask(SMOD_20M, SMOD_18M)
SMOD_21M<- mask(SMOD_21M, SMOD_18M)
NDVImax.2017 <- mask(NDVImax.2017, SMOD_17M)
NDVImax.2018 <- mask(NDVImax.2018, SMOD_17M)
NDVImax.2019 <- mask(NDVImax.2019, SMOD_17M)
NDVImax.2020 <- mask(NDVImax.2020, SMOD_17M)
NDVImax.2021 <- mask(NDVImax.2021, SMOD_17M)Nombre de pixels des couches SMOD:
## [1] 3501381 3501381 3501381 3501381 3501381
Nombre de pixels des couches NDVI:
## [1] 3501381 3501381 3501381 3501381 3501381
Le nombre de pixels est désormais identique pour l’ensemble des couches NDVI max et SMOD.
4.2.4 Réduction de la zone d’étude
Le nombre de pixels étant très important, le script requiert 6 heures de calcul pour être exécuté. Dans le cadre de notre projet étudiant, le choix a donc été de délimiter une zone d’étude restreinte afin de réduire le temps de calcul nécessaire.
4.2.4.1 Extraction d’une portion des jeux de données raster
Une zone a été définie manuellement pour couvrir à la fois la vallée et les versants et a été appliquée aux rasters NDVI et SMOD, afin d’extraire une portion de ces rasters selon un modèle d’étendue.
cropbox2 <-c(336000,339000,5083000,5086000)
SMODcrop17 <- crop(SMOD_17M, cropbox2)
SMODcrop18 <- crop(SMOD_18M, cropbox2)
SMODcrop19 <- crop(SMOD_19M, cropbox2)
SMODcrop20 <- crop(SMOD_20M, cropbox2)
SMODcrop21 <- crop(SMOD_21M, cropbox2)
NDVIcrop17 <- crop(NDVImax.2017, cropbox2)
NDVIcrop18 <- crop(NDVImax.2018, cropbox2)
NDVIcrop19 <- crop(NDVImax.2019, cropbox2)
NDVIcrop20 <- crop(NDVImax.2020, cropbox2)
NDVIcrop21 <- crop(NDVImax.2021, cropbox2)4.2.4.2 Visualisation du découpage appliqué
Figure 5: Visualisation de la couche SMOD 2017 après découpage selon une zone délimitée
4.2.5 Préparation du tableau de données NDVI max
Les données de NDVI max sont récupérées pour chaque année et pour chaque pixel.
NDVI_all <- data.frame(NDVI2017=NDVIcrop17[which(!is.na(NDVIcrop17[]))],NDVI2018=NDVIcrop18[which(!is.na(NDVIcrop18[]))],NDVI2019=NDVIcrop19[which(!is.na(NDVIcrop19[]))],NDVI2020=NDVIcrop20[which(!is.na(NDVIcrop20[]))],NDVI2021=NDVIcrop21[which(!is.na(NDVIcrop21[]))])
head(NDVI_all)| NDVI2017 | NDVI2018 | NDVI2019 | NDVI2020 | NDVI2021 |
|---|---|---|---|---|
| 0.3653489 | 0.3707226 | 0.3768148 | 0.4216587 | 0.4651088 |
| 0.3420864 | 0.3607130 | 0.3040986 | 0.3711887 | 0.4289954 |
| 0.1543539 | 0.2082854 | 0.0792078 | 0.1866200 | 0.1487487 |
| 0.0955285 | 0.0887893 | 0.0729308 | 0.1797512 | 0.0875441 |
| 0.2148243 | 0.2155243 | 0.2773324 | 0.2985634 | 0.2387886 |
| 0.3455674 | 0.3172919 | 0.4380009 | 0.4517799 | 0.4957235 |
4.2.6 Préparation du tableau de données SMOD
De même, les données de SMOD sont récupérées pour chaque année et pour chaque pixel.
SMOD_all <- data.frame(SMOD17M = SMODcrop17[which(!is.na(SMODcrop17[]))], SMOD18M = SMODcrop18[which(!is.na(SMODcrop18[]))], SMOD19M = SMODcrop19[which(!is.na(SMODcrop19[]))], SMOD20M = SMODcrop20[which(!is.na(SMODcrop20[]))], SMOD21M = SMODcrop21[which(!is.na(SMODcrop21[]))])
head(SMOD_all)| SMOD17M | SMOD18M | SMOD19M | SMOD20M | SMOD21M |
|---|---|---|---|---|
| 145.000 | 161.5 | 140.25 | 162.4375 | 116 |
| 145.000 | 165.5 | 143.25 | 163.7500 | 116 |
| 145.000 | 162.5 | 145.75 | 165.2500 | 116 |
| 149.125 | 161.0 | 149.50 | 166.0000 | 116 |
| 157.375 | 161.0 | 154.50 | 166.0000 | 116 |
| 157.375 | 160.5 | 153.25 | 163.5625 | 116 |
4.3 Calcul des pentes de tendance
La suite du projet consiste à calculer les pentes de tendance sur les 5 années à la fois pour les données NDVI max et SMOD. Le jeu de données est très faible (seulement 5 années) mais l’objectif de ce projet était surtout de définir une méthode d’analyse bien que celle ci soit plus adaptée à une série temporelle plus longue.
4.3.1 Calcul des pentes pour le NDVI max
4.3.1.1 Boucle de calcul des pentes et des pvalue associées
Greening <- data.frame(Pente = rep(NA, nrow(NDVI_all)),
Pvalue = rep(NA, nrow(NDVI_all)))
for(n in 1:nrow(NDVI_all)){
ndvi_pixel = data.frame(NDVI = as.numeric(NDVI_all[n,]), Year = c(2017:2021))
Pente <- summary(mblm(NDVI ~ Year, data=ndvi_pixel))$coefficients[2,1]
Pvalue <- summary(mblm(NDVI ~ Year, data=ndvi_pixel))$coefficients[2,4]
Greening[n,"Pente"] <- Pente
Greening[n,"Pvalue"] <- Pvalue
}4.3.1.2 Génération du raster des pentes NDVI max
NDVI_pente <- NDVIcrop17
NDVI_pente[which(!is.na(NDVI_pente[]))]=Greening[,"Pente"]Figure 6: Raster des pentes associées à l’évolution du NDVI max entre 2017 et 2021
Comme attendu, la plupart des pentes sont proches de 0 traduisant une faible évolution des données NDVI max sur les 5 années considérées.
4.3.2 Calcul des pentes pour le SMOD
4.3.2.1 Boucle de calcul des pentes et des pvalue associées
Melting <- data.frame(Pente = rep(NA, nrow(SMOD_all)),
Pvalue = rep(NA, nrow(SMOD_all)))
for(n in 1:nrow(SMOD_all)){
SMOD_pixel = data.frame(SMOD = as.numeric(SMOD_all[n,]), Year = c(2017:2021))
Pente <- summary(mblm(SMOD ~ Year, data=SMOD_pixel))$coefficients[2,1]
Pvalue <- summary(mblm(SMOD ~ Year, data=SMOD_pixel))$coefficients[2,4]
Melting[n,"Pente"] <- Pente
Melting[n,"Pvalue"] <- Pvalue
}4.3.2.2 Génération du raster des pentes SMOD
SMOD_pente <- SMODcrop17
SMOD_pente[which(!is.na(SMOD_pente[]))]=Melting[,"Pente"]Figure 7: Raster des pentes associées à l’évolution du SMOD entre 2017 et 2021
Les valeurs de pente positives sont plutôt liées aux zones de vallée alors que les valeurs négatives sont plutôt liées aux versants traduisant surtout une précocité de la fonte des neiges sur les versants.
4.4 Etude de la corrélation entre l’évolution du NDVI max et l’évolution du SMOD
Un modèle linéaire est appliqué aux deux variables quantitatives.
##
## Call:
## lm(formula = df$GreeningPente ~ df$MeltingPente)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.123012 -0.005490 -0.000536 0.005438 0.140579
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.442e-03 5.543e-05 44.06 <2e-16 ***
## df$MeltingPente -1.808e-04 4.764e-06 -37.95 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01646 on 89998 degrees of freedom
## Multiple R-squared: 0.01575, Adjusted R-squared: 0.01574
## F-statistic: 1440 on 1 and 89998 DF, p-value: < 2.2e-16
Figure 8: Relation entre NDVI max et SMOD
Une très faible relation négative attendue est observée entre l’évolution du SMOD et l’évolution du NDVI max.
Cependant, le modèle a été appliqué à l’ensemble des pixels sans tri préalable des pentes en fonction de la p value associée.
La relation est à nouveau testée seulement pour les pixels qui présentent les pvalue les plus faibles (p value < 0.07).
##
## Call:
## lm(formula = df2$GreeningPente ~ df2$MeltingPente)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.089626 -0.007452 -0.000791 0.008617 0.082001
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.308e-03 2.632e-04 16.37 <2e-16 ***
## df2$MeltingPente -3.759e-04 1.919e-05 -19.58 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01889 on 5802 degrees of freedom
## (1689 observations effacées parce que manquantes)
## Multiple R-squared: 0.062, Adjusted R-squared: 0.06184
## F-statistic: 383.5 on 1 and 5802 DF, p-value: < 2.2e-16
Figure 9: Relation entre NDVI max et SMOD pour les p values faibles
Une relation négative est observée entre l’évolution du SMOD et le
NDVI max (p value < 0.01). Ainsi, des valeurs de SMOD faibles (dates
de fonte précoces) seraient liées à des valeurs de NDVI élevées.
Ce
phénomène peut être expliqué par des périodes végétatives plus longues
permises par une fonte précoce de la neige qui favoriseraient le
verdissement de ces zones.
Cependant, le R2 associé à la relation
linéaire est très faible (0.06), probablement en raison du nombre de
données élevé issues de zones très diversifiées qui peut augmenter la
variabilité observée et du faible nombre d’années étudiées qui ne permet
pas d’avoir des évolutions significatives. En effet, peu de variations
du NDVI max étaient observées sur les 5 années considérées.
Il
faudrait appliquer ce modèle à une série temporelle plus longue ou à une
zone d’étude plus homogène (exposition, altitude, type d’habitats) pour
vérifier l’existence d’une corrélation entre les deux variables.